function gaussian_elimination()
% GAUSSIAN_ELIMINATION 高斯消元法演示
% 
% 演示高斯消元法的基本原理和实现

fprintf('=== 高斯消元法 ===\n\n');

%% 1. 基本高斯消元法
fprintf('1. 基本高斯消元法\n');
fprintf('求解线性方程组 A*x = b\n\n');

% 创建测试矩阵
A = [2, 1, -1; 
     -3, -1, 2; 
     -2, 1, 2];
b = [8; -11; -3];

fprintf('原始矩阵 A =\n');
disp(A);
fprintf('右端向量 b = [%g; %g; %g]\n', b);

% 执行高斯消元
[A_upper, b_upper, multipliers] = gauss_elimination_basic(A, b);

fprintf('消元后的上三角矩阵:\n');
disp(A_upper);
fprintf('消元后的右端向量: [%g; %g; %g]\n', b_upper);

% 后向替换求解
x = backward_substitution(A_upper, b_upper);
fprintf('解: x = [%.6f; %.6f; %.6f]\n', x);

% 验证解
residual = A * x - b;
fprintf('残差: ||A*x - b|| = %e\n', norm(residual));

%% 2. 高斯消元的矩阵形式
fprintf('\n2. 高斯消元的矩阵形式\n');
fprintf('高斯消元等价于 L*A = U，其中L是单位下三角矩阵\n');

% 构造消元矩阵
n = size(A, 1);
L_matrices = cell(n-1, 1);
L_inv_matrices = cell(n-1, 1);

A_temp = A;
for k = 1:n-1
    % 第k步消元矩阵
    L_k = eye(n);
    for i = k+1:n
        L_k(i, k) = -A_temp(i, k) / A_temp(k, k);
    end
    L_matrices{k} = L_k;
    L_inv_matrices{k} = inv(L_k);
    
    fprintf('第%d步消元矩阵 L_%d =\n', k, k);
    disp(L_k);
    
    % 应用消元矩阵
    A_temp = L_k * A_temp;
end

% 计算总的L矩阵
L_total = eye(n);
for k = 1:n-1
    L_total = L_total * L_inv_matrices{k};
end

fprintf('总的L矩阵 (L = L_1^(-1) * L_2^(-1) * ...) =\n');
disp(L_total);
fprintf('上三角矩阵 U =\n');
disp(A_temp);

% 验证 L*U = A
LU_product = L_total * A_temp;
fprintf('验证 L*U - A 的误差: %e\n', norm(LU_product - A, 'fro'));

%% 3. 高斯消元失败的情况
fprintf('\n3. 高斯消元失败的情况\n');
fprintf('当主元为零时，高斯消元会失败\n');

% 创建主元为零的矩阵
A_singular = [0, 1, 2; 
              1, 2, 3; 
              2, 3, 4];
b_singular = [3; 6; 9];

fprintf('奇异矩阵 A =\n');
disp(A_singular);
fprintf('矩阵的秩: %d\n', rank(A_singular));
fprintf('矩阵的行列式: %g\n', det(A_singular));

try
    [A_fail, b_fail] = gauss_elimination_basic(A_singular, b_singular);
    fprintf('消元成功完成\n');
catch ME
    fprintf('消元失败: %s\n', ME.message);
end

%% 4. 数值不稳定性示例
fprintf('\n4. 数值不稳定性示例\n');
fprintf('小主元导致的数值不稳定\n');

% 创建小主元矩阵
epsilon = 1e-10;
A_unstable = [epsilon, 1; 
              1, 1];
b_unstable = [1; 2];

fprintf('小主元矩阵 (ε = %e):\n', epsilon);
disp(A_unstable);

% 不选主元的消元
[A_unstable_elim, b_unstable_elim] = gauss_elimination_basic(A_unstable, b_unstable);
x_unstable = backward_substitution(A_unstable_elim, b_unstable_elim);

fprintf('不选主元的解: x = [%.6f; %.6f]\n', x_unstable);

% 精确解
x_exact = A_unstable \ b_unstable;
fprintf('精确解: x = [%.6f; %.6f]\n', x_exact);
fprintf('误差: ||x_computed - x_exact|| = %e\n', norm(x_unstable - x_exact));

%% 5. 运算量分析
fprintf('\n5. 运算量分析\n');
fprintf('高斯消元的运算量为 O(n^3)\n');

% 测试不同规模的运算时间
sizes = [50, 100, 200, 400];
times = zeros(size(sizes));
flops_theoretical = zeros(size(sizes));

fprintf('\n矩阵规模    实际时间    理论运算量\n');
fprintf('--------    --------    ----------\n');

for i = 1:length(sizes)
    n = sizes(i);
    A_test = randn(n, n) + eye(n);  % 确保非奇异
    b_test = randn(n, 1);
    
    tic;
    [A_elim, b_elim] = gauss_elimination_basic(A_test, b_test);
    x_test = backward_substitution(A_elim, b_elim);
    times(i) = toc;
    
    flops_theoretical(i) = 2*n^3/3 + n^2;  % 理论运算量
    
    fprintf('%6d      %6.3f s    %8.0f\n', n, times(i), flops_theoretical(i));
end

% 绘制时间复杂度图
figure('Name', '高斯消元时间复杂度');
loglog(sizes, times, 'bo-', 'LineWidth', 2, 'MarkerSize', 8);
hold on;
loglog(sizes, flops_theoretical/1e9, 'r--', 'LineWidth', 2);
xlabel('矩阵规模 n');
ylabel('时间 (秒) / 运算量 (10^9)');
title('高斯消元时间复杂度');
legend('实际时间', '理论运算量/10^9', 'Location', 'northwest');
grid on;

end

function [A_upper, b_upper, multipliers] = gauss_elimination_basic(A, b)
% GAUSS_ELIMINATION_BASIC 基本高斯消元法
% 
% 输入:
%   A - n×n 系数矩阵
%   b - n×1 右端向量
% 输出:
%   A_upper - 上三角矩阵
%   b_upper - 消元后的右端向量
%   multipliers - 消元乘数

n = size(A, 1);
A_upper = A;
b_upper = b;
multipliers = zeros(n, n);

for k = 1:n-1
    % 检查主元是否为零
    if abs(A_upper(k, k)) < eps
        error('主元为零，需要选主元策略');
    end
    
    % 消元
    for i = k+1:n
        multiplier = A_upper(i, k) / A_upper(k, k);
        multipliers(i, k) = multiplier;
        
        % 更新矩阵行
        A_upper(i, k:n) = A_upper(i, k:n) - multiplier * A_upper(k, k:n);
        b_upper(i) = b_upper(i) - multiplier * b_upper(k);
    end
end

end

function x = backward_substitution(U, b)
% BACKWARD_SUBSTITUTION 后向替换算法
% 求解上三角系统 U*x = b

n = size(U, 1);
x = zeros(n, 1);

for i = n:-1:1
    sum_term = 0;
    for j = i+1:n
        sum_term = sum_term + U(i, j) * x(j);
    end
    x(i) = (b(i) - sum_term) / U(i, i);
end

end